SPLAY STATES IN FINITE PULSE-COUPLED NETWORKS OF 
EXCITABLE NEURONS 
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Abstract. The emergence and stability of splay states is studied in fully coupled finite networks 
of N excitable quadratic integrato-and-fire neurons, connected via synapses modeled as pulses of 
finite amplitude and duration. For such synapses, by introducing two distinct types of synaptic 
events (pulse emission and termination), we were able to write down an exact event-driven map for 
the system and to evaluate the splay state solutions. For M overlapping post synaptic potentials 
the linear stability analysis of the splay state should take in account, besides the actual values of the 
membrane potentials, also the firing times associated to the M previous pulse emissions. As a matter 
of fact, it was possible, by introducing M complementary variables, to rephrase the evolution of the 
network as an event-driven map and to derive an analytic expression for the Floquet spectrum. We 
find that, independently of M, the splay state is marginally stable with N — 2 neutral directions. 
Furthermore, we have identified a family of periodic solutions surrounding the splay state and sharing 
the same neutral stability directions. In the limit of 5-pulses, it is still possible to derive an event- 
driven formulation for the dynamics, however the number of neutrally stable directions, associated 
to the splay state, becomes A'. Finally, we prove a link between the results for our system and a 
previous theory [Watanabe and Strogatz, Physica D, 74 (1994), PP- 197- 253] developed for networks 
of phase oscillators with sinusoidal coupling. 
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1. Introduction. The dynamics of networks made up of many elements with a 
high degree of connectivity is often studied in the infinite size Hmit, which allows to 
apply approaches borrowed from statistical physics. In particular, for globally coupled 
neural networks this amounts to finding the distribution of the membrane potentials 
satisfying a Fokker-Planck equation with specific boundary conditions corresponding 
to the spike emission and reset of the neurons [H [T]. In contrast, general techniques 
to deal with the dynamics of finite size ensembles are not yet fully developed, not 
even for the analysis of the linear stability of periodic solutions. 

In this paper we investigate the stability of splay states (also known as antiphase 
states or "ponies on a merry-go-round" ) [TSl [5] ■ In a splay state all the N elements 
follow the same periodic dynamics x{t) {x{t + N ■ T) = x{t)) but with different time 
shifts evenly distributed at regular intervals AT — kT, with k — 1, . . . ,N. Exper- 
imental observations of splay states have been reported in multimode laser systems 
[30j and electronic circuits [-f!. Numerical and theoretical analyses have been de- 
voted to splay states in Josephson junction arrays [121 [551 UHl 13], globally coupled 
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Ginzburg-Landau equations [16], globally coupled laser models [23], trafSc models 
[21], and pulse-coupled neuronal networks [T]. In the latter context, splay states have 
been usually investigated for leaky-integrate-and-fire (LIF) neurons and in general for 
neuronal models which can be assimilated to phase oscillators (rotators) [Tl 
The first detailed stability analysis of LIF neuron oscillators was performed by devel- 
oping a mean- field approach in the infinite network limit [U [28] . Finite size stability 
analysis for supra-threshold neurons, namely for LIF |32j and for generalized neuronal 
models in the oscillatory regime [B], have been more recently developed, based on the 
linearization of a suitable Poincare map. 

The model analyzed in this paper is a a fully coupled network of excitable neurons, 
governed by the quadratic integrate- and-fire equation (QIF). The QIF is the canonical 
model for type I neuronal excitability as it is the quadratic normal form for the saddle- 
node invariant cycle (SNIC) bifurcation [9 . The neurons are coupled with positive 
pulses, modeling excitatory synapses. We focus our analysis on the persistent activity 
of the network that is induced by the recurrent excitation and that co-exists with an 
inactive ground state. 

Analyzing this type of activity is of significant relevance to neuroscience. The 
bistable sustained activity has been shown to be the neuronal basis for working mem- 
ory [TU] [H] . Such self-sustained elevated activity has been recorded in delayed re- 
sponse tasks where the memory trace must be retained in order to generate appropri- 
ate responses. Furthermore, so-called cortical up-states, observed during anesthesia 
and during sleep, are also considered to be generated by the intrinsic excitatory synap- 
tic connectivity with the constituent neurons being excitable (as opposed to intrinsic 
oscillators) . There are also indications that these sustained up-states are largely asyn- 
chronous. In fact, theoretical studies have suggested that asynchrony is a requirement 
for stable maintenance of synaptically sustained neural activity [20] • 

Furthermore, previous computational work proposed that perturbing the asyn- 
chronous structure of the sustained activity leads to its destabilization [TU [5] . It is 
thus important to determine specifically the stability and the structure of the asyn- 
chronous sustained activity. This item has been addressed in the infinite size limit 
within the mean-field approximation [TH| [T7] and the the role of asynchrony and syn- 
chrony in sustained neural activity has been studied for a pair of neurons [14) . However 
sustained cortical activity appears to be generated by local circuits in the cortex, i.e. 
networks with a limited number of neurons. Hence in our work we seek to understand 
the stability of asynchronous activity self-sustained by a finite size network. 

In this paper, as already mentioned, we analyze the splay states, which represent 
highly symmetric states. We perform an analytical linear stability analysis of the splay 
states for finite size networks when the post-synaptic potentials (PSPs) are modeled 
as square pulses of finite amplitude and duration. We focus on fast excitatory synaptic 
coupling as a basic mechanism to generate the reverberative self-sustained activity. 
This corresponds to AMPA receptor-mediated glutamatergic synapses that have a 
typical decay-time constant of about 5 msecs [25] . Traditionally such synapses are 
modeled as a double exponential function (or an a-function) with a finite rise-time 
and a decay-time governed by the synaptic time-constant [12 . 

Here we use a simpler version of this model: we keep the idea of the characteristic 
synaptic time scale while leaving aside the dynamics by modeling the synaptic currents 
as square pulse steps. The advantage of such a minimal model is that it makes the 
network dynamics tractable for our analysis, while giving us control over the synaptic 
duration. 
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In order to study the finite size network, we derive an event driven map for 
the evolution of the membrane potentials of the neurons, by introducing two kinds 
of synaptic events: synaptic pulse emission and termination. This approach allows 
us to derive an analytic, but implicit, expression for the splay state for two kinds of 
synaptic models: step and (5-pulses. Furthermore, the linear stability analysis requires 
the investigation of the linearized dynamics of the model. It should be mentioned 
that memory effects should be taken in account whenever the duration of the post- 
synaptic potentials last sufficiently to lead to overlaps among the emitted pulses. For 
M overlapping pulses, the linearized dynamics can be rewritten as an event-driven 
map by including Af additional variables. This at variance with the usual approach, 
where the memory effect due to the linear super-position of a- or exponential pulses 
emitted in the past is taken in account by a self-consistent field [Tl. Finally, by 
employing the event-driven formulation we have analytically obtained the Floquet 
spectra associated to the splay state for step and 5-pulses. 

The paper is organized as follows. In §2 the model and the possible dynamical 
regimes are introduced. The event-driven map for step and 5-pulses is derived in §3, 
while the linear stability analysis of splay states is performed in §4 for step pulses and 
in §5 for (5-pulses. §6 is devoted to the description of other periodic states observable in 
the present model. Finally, in §7 the results are summarized and discussed. Analytical 
expression for the firing rates of the splay states in small networks are reported in §A. 
Furthermore, in §B we report an analytical expression for the splay state membrane 
potentials derived in the continuum limit. §C contains a formal proof for our model, 
in the case of non-overlapping pulses, that the Floquet spectrum associated with the 
splay states contains N — 2 marginally stable directions. 

2. Model and Dynamical Regimes. In this Section, we will introduce our 
model and the specific state which is the main subject of investigation of our analysis, 
namely the splay state. In particular, we consider a pulse-coupled fully connected 
excitatory network made of quadratic integrate and fire neurons, whose dynamics is 
governed by the following equation 



where the n-th spike is emitted at time i„, once the neuron reaches the threshold 
value Vi(t~) = oo; afterwards it is immediately reset to the value Vi(t^) = —oo. For 
a constant synaptic current / < 1, the neuron has a stable fixed point at Vrest = 
— -y/r^-^ and an unstable ones at Vu = +^/^ — I- The dynamics is excitable with w„ 
representing the threshold to overcome to observe "an excursion" towards infinity (a 
spike) before relaxing to the rest state at Vrest • This amounts to saying that if the 
initial value Vi(t = 0) < Vrest also at all the successive times the membrane potential 
will remain smaller than Vrest- While if Vrest < Vi{t = 0) < w„ the membrane potential 
will tend asymptotically to Vrest- Furthermore, for / > 1 the neuron fires periodically 
with frequency = \jl — l/(7rT). 

Since the network is fully connected, with equal synaptic weights, all neurons 
receive the same synaptic current l(t) that is the linear superposition of all the pulses 
emitted in the network up to the time t. In particular, as schema for the Post-Synaptic 
Potentials (PSPs) we consider step functions of finite duration Tg and amplitude 
J = G/{NTs), therefore the current reads as: 



r 



dt 



v1-l + I{t) z = 



N 



(2.1) 




(2.2) 



{tn} 
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where Q(x) is the Heaviside function, the sum runs over ah the spike times t„ < t, 
and the coupUng is normahzed by the number of neurons N to ensure that the total 
synaptic input will remain finite in the limit N ^ oo. We have consider pulses of the 



form (2.2 1 as the simplest example of PSPs allowing to take in account spatial and 
temporal summation of stimuli, due to their finite duration and amplitude. 

In the limit Tg — >■ 0, the PSPs will become (5-pulses and in this case the synaptic 
current can be rewritten as follows: 

m = ^J2^(t'*"^ (2.3) 

By following |18| . we can derive the average firing rate v in infinite size network, 
in this case the spiking frequency of the single neuron is simply given by 



^ = (2.4) 

TTT 

where Gv is the total synaptic current received by each single neuron, this result is 



valid both for the step PSPs (2.2 1 as well as for the (5-pulses. By solving the implicit 
equation above one gets 



-1,2 = (2.5) 

therefore there are two branches of solutions, we will re-examine this point later. 
Let us just mention that these solutions have been associated with the asynchronous 
persistent states emerging in networks composed by inhibitory and excitatory QIF 
populations (18j . 

A particular example of asynchronous state emerging in globally coupled networks 
is the so-called splay state |5Hl El] ■ This regime is characterized by a sequential firing 
of all the neurons with a constant network interspike interval (NISI) T, while the 
dynamics of each neuron is periodic with period N ■ T. This state has been classified 
as an asynchronous regular state [S]. Splay states have been found in an all-to-all 
pulse coupled excitatory network for LIP models |32j as well as for general neuronal 
models [HI , and in inhibitory networks for (5-pulses [1] . 

3. Event-driven map. As previously done in P^[5^ for LIF neuronal models, 
we would like to derive an event-driven map for the setup considered in the present 
paper. The event-driven map gives the exact evolution of the system, described by 



the set of N ODEs (2.1) plus the variable describing the synaptic current, from an 
event to the successive one. Therefore the continuous time evolution is substituted 
by a map with discrete time. 



Let us first consider PSPs that are step pulses of duration Ts as reported in (2.2). 
In the last part of the section we will derive the event-driven map also in the (5-pulses 
limiting case. 

3.1. Step pulses. In the case of step pulses, two type of events should be dis- 
tinguished: pulse emission (PE) and pulse termination (PT). Both events induce 
an instantaneous change of the synaptic current by a constant value: the current 
will increase (resp. decrease) by a quantity J for PE (resp. PT). In order to inte- 
grate the system it is not sufficient to know the value of the membrane potentials 
and of the synaptic current at a certain time t. The system evolution will depend 
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also on the termination times of the previous pulses received by the neuron that are 
"active" (still contributing to the synaptic current) at time t. Therefore one needs 
to know the ordered list of the future PT times {Sj{t)} with j = 1,. . . ,K, where 
t < Si(t) < S2{t) < ■ ■ ■ < SK{t)- The number K{t) of these events is in general not 
constant and it represents the number of overlapping pulses at time t, which amounts 
to a synaptic current I{t) = K{t)J . Let us now discuss separately how the PE and 
PT events influence the neural dynamics in order to derive an event-driven map. 

Pulse Emission. Suppose that at time i„ the neuron q emits a spike and that at 
time t~ there were K overlapping pulses. One can obtain the value of the membrane 
potential for the neuron i at the next event, occurring at t„ + At, by integrating 



equation ([21]) with I{t) = {K + 1) J 

r..(t.+At) /•*"+^* dt 

dtt) X^ + {K+l)J-l ' 



(3.1) 



How to determine the time interval At will be explained in the following. Due to the 



simple form of the PSP we can integrate (3.1 1 analytically, obtaining 

r H{v,{t+),K+l,M), i^q 
w,(t„ + At)=<^ (3.2) 
[ H*{K+l,M), i = q 

with 

i/(x,X,t) = I^^^^«±^, H*{K,t)^-l/PK{t), (3.3) 
and with the function Pk defined as follows 



PK{t) - <^ 



KJ < 1 



KJ > 1 



tanh(Vl--fS'Jt/-r) 
Vl-KJ 



tan(V-fS'J-lt/-r) 
VK.J-1 



(3.4) 



Furthermore, the list of the future PT times should be updated by adding SK+i{tn) = 

Pulse Termination. Let us now consider a PT occurring at time tpx when there 
were K > I overlapping pulses present in the network. The membrane potential of 
the i-th neuron at the next event, occurring at tpr + At, can be obtained by solving 
the following integral 



X^ + {K-1)J-1 7,+ 



(3.5) 



which gives 



v,{tpT + At)^H{v,{t+j.),K~l,At) . (3.6) 

At each pulse termination the list of the PT times {Sj{tpT)} should be updated by 
throwing away the smallest time 5*1 and by relabeling the other times as S'j(tp^) = 
Sj+i{tprp) with j ~ 1, K — I. 



^Please notice that in the excitable case {KJ < 1) one gets a single valued function from the 
integral ( |3.1| due to the fact that, depending on the initial value of the membrane potential, the 
dynamics remains segregated in one of the three intervals Vi{t) < Vrest or Vrest < ^iit) < Vu or 
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Determination of the Integration Time-Lapse. After each event PE or PT at time 
t* , one should determine the time interval At until the next event. In particular, one 
should understand if the next event will be a PE or a PT. In order to resolve this 
dilemma, the next presumed firing time E{t*) occurring in the network has to be 
firstly determined on the basis of the values of the membrane potentials and of the 
synaptic current at time t* . In the absence of any intermediate event, since we are 
considering a fully coupled system, the neuron p with highest membrane potential 
value Vp{t*) is going to fire at time E{t*). This time can be determined by imposing 



that H{vp{t*),K, E{t*) ~ t*) ^ oo, with H given by eq. (3.3), namely 



KJ <l t* 



tanh 



-1 f Vl-Kj \ 



YJ 



E{t*) = < 



(3.7) 



KJ >l t* + 



VKJ~1 



tan-if^^ 



where K is the number of overlapping pulses immediately after the event at t* . In 
order to understand the type of the next event E{t*) should be compared with Si{t*) 
to determine which is the smaller one. li K = then At = E{t*) — t* automatically, 
otherwise 



At = mill {E{t*), Slit*)} -t* . (3.8) 
The event-driven map will be therefore a combination of the two above described 



integration steps. After each event the potential will be given by eq. (3.2) or eq. (3.61 
depending if the event is a PE or a PT. 

Co-moving frame. A further simplification to the above scheme can be obtained 
by exploiting the fact that for globally coupled networks the neuron firing order is 
preserved. Since the firing order is directly related to the membrane potential value, 
we can order sequentially the membrane potentials, i.e. vi{t) > V2{t) > ■ ■ ■ > VN{t), 
and introduce a co-moving frame. This amounts to relabeling the neuron closest-to- 
threshold as 1, and when it fires at time t„ to reset the potential value as vi{t^) 
z)Ar(t+) = — oo and to shift the indexes of all the others ? — >■ {i — 1) for i > 2. 



Furthermore, due to the reference frame transformation, Eq. (3.2 ) has to be modified: 
namely the evolution map should be rewritten as Vi{tn-\-At) — H{vi+i{t^), K-\-l, At), 
for i = 1, . . . , AT - 1 and t;Ar(t„ + At) = H*{K + l,At). 

3.1.1. Splay-state. For the splay state regime, the event-driven map outlined 
above simplifies noticeably and furthermore it can be explicitly written. The splay 
state is characterized by a constant NISI: T. Furthermore, due to the regular spike 
emission the PT times can be all written in function of Si (t) as Sj {t) = Si (t) + {j — 
1) ■ T. In general, it is useful to rewrite Tg as a function of T, as follows 

Ts = MT + To (3.9) 

where K ~ M is the number of overlapping PSPs just before the spike emission. 
To < T and let us define Ti = T — Tq. Please notice that for a splay state K can 
assume only two values, namely M and M -I- 1 as shown in Fig. |3.1| In the case of 
non-overlapping pulses: M = 0, Tg = Tq, and Ti =T — Ts. This case is illustrated in 



Fig. 3.1 
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In order to determine the value of the coupling Gm required to have exactly M 
overlapping pulses, let us employ, as a first approximation, the mean field equation 



MT, 



(2.5) with the condition vi = M/{NTs), which is equivalent to assuming that = 

(3.10) 



M 



G 



M 



NTs 2r27r2 
then we can invert the above equation and obtain the critical coupling 



Jm 



Gm 



I 



(3.11) 



If J < Ji there is no overlap among two successive emitted PSPs. When Jm < J < 
Jm+1, M pulses overlap. The synaptic current can take only the following two values 



lit) 




tn<t< tn + To 



(3.12) 



as clearly illustrated in Fig. |3.1[ In particular, if Tq = 0, we will have always exactly 
M overlapping pulses, since each PE will coincide with a PT, and / — MJ. 



(a), 



M = 

I 



T (n 
I 



2J 

J 



M = l 



I 



To'" 



IT 



t;- 



(c) 

I 

3J. 

2J ■ 
J ■ 



Sl{tn)t"+1 

M = 2 



tn+2 



I 



tn-1 



Sl(tnf"+^ Si (tn) tn+2 



I 



Tr 



tn Sl(tn) 



4. 



t t 

S2(tn) 



SsCtn) 



Fig. 3.1. PSPs in a splay state can overlap M times, (a) PSPs overlapping M = times (no 
overlaps); (b) PSPs overlapping M = 1 times; (c) PSPs overlapping M = 2 times. Independently 
of the value of M the synaptic current can take only two values in the time interval between two 
spikes: namely, during I{tn < t < tn + Tq"') = (M + 1) J and I(tn + Tq'"' <t < = MJ. 
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For the splay state we can rewrite the dynamics of each neuron i between two 
successive spikes occurring at t„ and as an exact map made of the following three 
steps. 

1. The first step starts with a PE at time t„, one can easily estimate the evolution 
of the membrane potential from time t+ to Ti when a PT will occur. Let us 
first define a;|"^ = Vi(t^) and — Vi(t^ + To), and order the membrane 
potentials as follows 



„(") 



in) 



> 



^ (») 

> X}^' = — oo 



(3.13) 



the last equivalence stems from the fact that a neuron has just fired and it 
has been reset. By employing the expression (3.2) one gets the following map 

(3.14) 



Fi(x|"\ro) - H{x'f'\M + l,To), I + N 
F^iTo)=H*{M + l,To), i^N 



(n) 



2/. 



with H and H* defined in ( |3.3| ). 
2. The second step corresponds to the integration of the equation of motion 
from the PT occurring at i„ + Tq and the time t~j^-^ immediately preceding 

the n + 1-th spike emission. By defining zf'^ — Wi(i^+i) and by employing 
equation (3.6) one gets 



An) 



(3.15) 



with H defined in (3.3). Due to the previous ordering, the next firing neuron 

in) 

will have the label 1, therefore z\ — oo and thus the denominator of the 
right side equation (3.15) should be zero: 

1 



By inserting (3.16) in (3.15) one gets: 



/ (n) {n)\ 



(AfJ-l)+2/(")y(") 



(«) 



(3.16) 



(3.17) 



3. The last step amounts simply to calculating the membrane potential change 
in going from t^^^^ to t^_^_i and introducing a co-moving frame to maintain 
the order among the membrane potentials also after each firing event. This 
amounts to writing 



P (An) \ _ Jn) 



for l<i< N ~1 



(3.18) 



and setting x^^~^^^ 



-oo. Since the event-driven map approach corresponds 
to a suitable Poincare section, we are left with — 1 variables, dropping the 
variable i = N. 

We can compute the complete event-driven map from spike time i„ to spike time 



tn+i, by combining the three above equations (3.14), (3.15) and (3.18) 



, (") 



, (") 
02 + a3xl^\ 



for 1 < i < iV - 1 



(3.19) 
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where the coefficients entering in ( |3.19| ) reads as: 



ao = (MJ - 1)/3m(Ti) + [{M + l)J - 1]/3m+i(To) 

ai - i - (MJ - l)(3M+iiTo)l3MiTi) 

a2 = l- [{M + 1) J - 1]I3m+i{To)Pm{Ti) 

an = -/3a/+i(T'o) - Pm{Ti) 



(3.20) 



Exact firing rate value.. In order to obtain the membrane potential values asso- 
ciated with the splay state one should impose that the splay state represents a fixed 
point for the event-driven map in the comoving frame, namely 



(3.21) 



Furthermore, once fixed G, TV, and Tg, one can determine the NISI T by solving 
iteratively equation ( 3.7 ) together with the set of equations for the membrane potential 
(3.19), with the requirement that x* = {xn = — oo) = -|-oo. Numerically, as a first 



guess for T we usually employ the mean- field result 1 /vi , given by the larger solution 
of (2.5 1. Then we evaluate the splay state by employing a bisection method to find 
the exact NISI. We stop the procedure whenever x* > 10^, with the constraint that 



the order (3.13) is maintained 



For a given set of parameters G, Tg and N we found at maximum two coexisting 
splay states (in agreement with the mean- field results). Beyond a minimal value of J, 
there is always one marginally stable splay state. When the splay states are two we 
found that the other one is unstable, as illustrated in the following in Fig. |4.4[ Let us 
stress that unstable branches of solutions exist only for non overlapping pulses (i.e. 
M = 0) as shown in Fig. 3.2 These numerical results will be confirmed by analytical 
analysis in §A for iV = 2, 3, 4 and J < Ji (M = 0). 

Notice that for iV = 2 only the marginally stable branch exists and the minimal 
firing rate reaches the value i/ = 0. Instead for > 2 the minimal firing rate of the 
marginally stable branch is ^ 0. The firing rates associated to the unstable branch, 
for finite N, reaches always the value = for some finite pulse amplitude J = J* . 
Finally, for TV — >■ cx) we have that J* (i^ = 0) -> oo. 

3.2. (5-pulses. In the case of (5-pulses, the formulation of the event-driven map is 
extremely simplified, since now there are only PE events. At the arrival of a (S-pulse, 
we can integrate eq. (2.1) with the current given by (2.3) between time t~ and t^, 
obtaining 



„(") 



yl ' = xl ' + Js for 1 < i < A - 1 , (3.22) 
G/{Nt). The evolution of the membrane potential in the time interval 



where J5 

t^ and i.^^x t)e easily obtained since it corresponds to eq. (3.6) with M 
- t^ = T, namely 



Ti = t; 



An) 



and 



(3.23) 



for i = 1, . . . , A^ — 1. Then we can combine eq. (3.22) and (3.23) with the change of 



reference frame (3.181 to obtain the corresponding event-driven map. The resulting 



map is identical to that found for the step function (3.19), apart the value of the 

-MT) + Js 



coefficients (3.20) that now become: 



ao 
ai 

0.2 

as 



1 

l-l3oiT)Js 
-MT) 



(3.24) 




Fig. 3.2. Frequencies of the splay states as a function of the synaptic strength J and for pulse 
duration TsN = 12 ms with r = 20 ms. Red line: N = 2, magenta line: N = 3, blue line: TV = 4, 
green line: N = oo. Black dotted lines separate regions with different number M of overlapping 
PSPs, Solid lines refer to the upper stable branches of the splay state. Dashed lines refers to the 
lower unstable branches of the splay state. For N = 2 the lower branch does not exist. 



Once fixed Js and TV and Ts, similarly to the case of step pulses, one can determine 
T together with the membrane potential values associated with the splay state by 



solving iteratively equation (3.7), and by applying iteratively the map (3.191 with 



—oo. The solution is numerically achieved 
oo (namely, x* > 10®) and the condition (3.131 is 



coefficients (3.241 starting from 
whenever x* — {xj^i ~ —oo) 
satisfied. 

We want to conclude this Section by mentioning the fact that in the limit N ^ oo 
we were able to derive an explicit analytic expression for the membrane potentials 
corresponding to a splay state. The detailed calculations are reported in §B. 

4. Linear Stability Analysis for Step Pulses. We are interested in the linear 
stability of the splay state in the case of step pulses for finite system size N. It is 
therefore useful to introduce the following vector notation for the membrane potentials 
at spike time 



n)\ 



(4.1) 



Furthermore, if we have more than one overlapping pulse, i.e. if M > 0, the actual 
state of the network will depend not only on the membrane potential values but also 
on the past M spike times {tk} with k = n — M, n — M + 1, . . . ,n — 1. However the 
formulation of the tangent space dynamics can be made simpler by introducing the 
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related time intervals r = t„ — t„ 



(n) (n) 



(«) 



(4.2) 



In this notation the splay state is a fixed point of the network dynamics satisfying 
the following relationships: 



X = Fsii) = F3{F2{y)) - F3{F2{F,{Sc))) , 



and 



Tj =j-T j = 1, 



(4.3) 



(4.4) 



4.1. Linearized Poincare map. In order to derive the equations of evolution 
in the tangent space for our case it is convenient to consider separately the three steps 
in eq. (4.3), please notice that now Tq"^ and T^"^ depend on the spike sequence index 



n since, for the perturbed dynamics these quantities are no more constant. 



Let us start by perturbing eq. (3.14): 



(4.5) 



with STn^^ = if Af = 0; where the coefficients are: 



dx 



(n) 



(n) 



l + p/ + l)J-l]/3i,+i(fo) 



(l-/3M+l(ro)iO' 



{M + l)J -l + ij 

(l-/3M+i(To)i^)2cos2 



(4.6) 



1 1 

V(m + l)g - ifo/r) ^ 



(4.7) 



sn 



Pli+iin) cos2 (V(M + l)J-lfo/r) ^ 



(4.8) 



As a second step we perturb F2 given by eq. (3.17), obtaining 



1 = 1,. ..,N 



(4.9) 



remember that if M = then Sy'"^'' = 0. The coefficient hi and ki are defined as: 
9i^2(2/i"\j/f^) 

yi,yi 



h, = 



dy\ 



(n) 



MJ -1 + yf 

{yi -yiY 



(4.10) 



hi — 



dF2{yr.y. 



(") 



yt 



dy, 



in) 



yi,Vi 



{yi - mf 



(4.11) 



12 



DIPOPPA, KRUPA, TORCINL AND GUTKIN 



Finally, the linearized equations associated with the reference frame change can 

(4.12) 



be obtained by perturbing eq. (3.181: 



,N-l 



please notice that Sx^^^ = due to the fact that in the comoving frame x 



(n) ^ 
N — 



therefore the evolution in the tangent space should deal with only — 1 perturbations 
associated to the membrane potentials. 

in] 

Then we need to compute how the time interval Tq is modified by the pertur- 

(n) 

bations, when M > 0. The key point here is that Tq depends on the previous spike 
times as follows: 



rp(n) ^ rp _ If _ f 



M. 



in) 
M 



(4.13) 



M ■ 



apparently one could be lead to think that we need only an extra variable: r 

(n) 

However, rjjj depends on all the M previous spike times and therefore we need to 

take in account also the perturbations of the other M—1 variables, namely rj"^^ a/- i • 
To obtain the evolution equations for these auxiliary M variables, let us consider 
the following relations 



_{n+l) 



in) 



in) 



and 



Jn+l) 



(n+1) 



(n) _ rp(n) 
^3-1 - ^0 



-n 



(4.14) 



(4.15) 



From (4.131 wc obtain the relation 5t 



equations (4.14 



in) 
M 



-st; 



in) 



From this relation and from 



and (4.15) (for positive M) we can easily obtain the evolution maps 



for the perturbed quantities: 



(") 



5t. 



(n+l) 
J=2,. 



.M 



6Tl 



(") 



St 



St. 



in) 
M 

(») 

i-1 



St 



(n) 
M 



(4.16) 



We are left just with the determination of St["'^ , this can be derived by remember 
ing that the time from the last FT until the next PE can be calculated by employing 
eq. (I3J| with K = M, Vp{t+ - n), and E(t*) - t* ^ 



in) 



G{y\ 



MJ < 1 



MJ > 1 



VMJ-1 



tanh 



tan 



n/MJ-1 



(4.17) 



and 



and we can obtain: 



dG 



w = 



dy\ 



(n) 



\MJ -l\ 



St[^^ = w5y\ 



in) 



(4.18) 



(4.19) 
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By combining eq. (4.5), (4.9 1, (4.12 1, (4.16), and (4.19), the complete map evo 



lution in the tangent space can be finahy written as follows: 



Sx 



Sx 



(n+l) 

i=l,...,N-2 

(n+1) 
N-1 



= Pi+iSx]^"^' + qi+iSxll{ + Ui+iSr 



M 



(4.20) 



— riSxl" + r2ST 



M 



St, 



(n+l) 



J=2,...,M 



{hiSi + kiSi), Ti = wdi, and r2 



where we have set pi = hidi, qi = k^di, Ui - 

-(l+WSl). 

In order to determine the stability of the splay state we should compute the 
Floquet spectrum by setting 



Sx 



Sri 



(n+l) 
N-1 
(n+l) 



V 5' 



M 



( sx^r \ 



"•^N-i 

X (") 
St{ ' 



\ 



(4.21) 



where /i; — e''''+"^' {I ~ 1, . . . , N + M — 1) are the so called (complex) Floquet 
multipliers, while A; (resp. uji) are real numbers termed Floquet exponents (resp. 
frequencies). If ||^;|| < 1 VZ (resp. ||^fc|| > 1 for at least one k) the splay state is 
stable (resp. unstable). Whenever the largest modulus of the Floquet multipliers is 
exactly one the system is marginally stable. 

The Floquet spectrum can be obtained by solving the following characteristic 



polynomial, obtained from eq. (4.20) 



P2^^l 



,N-2 



M 



,N-2 



J2k=3 "fc 



9j Mi 



.JV-fc 



(4.22) 



which admits N + M ~ 1 solutions. 



4.2. Floquet multipliers. As stated by Watanabe and Strogatz (WS) 29] for a 
network on N fully coupled phase oscillators with sinusoidal coupling, the system has 
in general iV — 3 marginally stable directions, furthermore for a splay state, which is 
a periodic solution, these directions reduce to N — 2. Therefore since also our model, 
as detailed in §C, satisfies the hypothesis for which the WS results apply, and since 
in the event-driven map formulation one degree of freedom is lost, we expect that for 
the splay states at least iV — 3 Floquet multipliers lie on the unit circle, as shown in 
Fig. |4.1| for M = 0. Furthermore, in presence of overlaps, i.e. for M > 0, the Floquet 
exponents associated to the auxiliary variables r'"-' do not influence the stability of 
the splay state, since these additional M exponents are located within the unit circle, 



and therefore associated to stable directions as shown in Fig. 4.2 and 4.3 
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-1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 

Re(n) Re((i) Re(n) 

Fig. 4.1. Floquet multipliers {fii} for the case with no overlap, i.e M = 0; (a) N = 3, 
marginally stable eigenvalue; (b) N = i, 1 marginally stable eigenvalue; (c) N = 8, 5 marginally 
stable eigenvalues. In this case we fixed J = 15 and NTs = 16 ms and we vary the network size. 



(a)l.o 




Fig. 4.2. Floquet multipliers {ni} for overlapping pulses, i.e. M > : (a) J = 15, M = 0, 2 
neutrally stable eigenvalues; (b) J = 25, M = 1, 2 neutrally stable eigenvalues; (c) J = 100, M = 6, 
2 neutrally stable eigenvalues. We have considered N = 5 and Ts = 3.2 ms. 



It is interesting to notice how the additional exponents associated to the auxihary 
variables emerge by increasing the number of overlaps. In particular, the number of 
overlaps can be increased from M to M + 1 by varying the coupling J from below to 
above the threshold Jm+i- At the threshold Jm+i a new variable tm+i is added to the 
event-driven map describing the system. Therefore the Floquet spectrum associated 
with the corresponding splay state solution has one additional eigenvalue. This new 
direction emerges as superstable at J = Jm+i being associated to a zero Floquet 



multiplier, as shown in Fig. 4.3 By further increasing J the new eigenvalue increases 
its modulus, which however remains always smaller than one. 



In Fig. 4.4 we report the Floquet multipliers associated to the unstable branch 
of splay state solutions, which coexist with the marginally stable branch for N > 2, 
as already mentioned in Sect. Ill A. 

5. Linear Stability for (5-pulses. In the case of (5-pulses the stability of the 
splay state can be inferred by theoretical arguments based on the symmetry of the 
considered model and of the specific pulse coupling. It is evident that the QIF 
model (2.1) for time symmetric pulses has a time reversal symmetry. This can 
be appreciated as follows. Given a solution v(t) — {vi{t), . . . ,VN{t)} we define 
w(t) — {wi{t), . . . ,wpf{t)} — —{vN{—t),...,vi{—t)}. It is clear from the time re- 



versal property of (2.1 1 that w(t) is a solution in between two spike emissions. Let us 



analyze if the symmetry is maintained also during spike emission, in the usual case 
vi will reach oo, then it will be reset to —oo and a constant value Js will be added 
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Fig. 4.3. (a) Floquet spectrum of the splay state in the complex plane for TsN = 16 ms, 
N = 5, in this case ,}\ = 16.42. Blue stars correspond to M = when J = 10.42 < Ji (blue), and 
J = 14.42 < Ji (red); and to M = 1 when J = 18.42 > Ji (cyan), and J = 22.42 > Ji (magenta) 




-1.0 -0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 1.5 2.0 -1 1 2 

Re((i) Re(!x) Re((,i) 



Fig. 4.4. Floquet spectrum of the splay state in the complex plane for the unstable branch, 
TsN = 16 ms, N = 10; (a) J = 8, (b) J = 10, (c) J = 12. 



to all the other membrane potentials. The membrane potential wi(t) reaching oo is 
equivalent to ujv(— reaching — oo. Backwards in time the reset and coupling consists 
of setting w^v to oo and subtracting Js from the other variables. Due to the minus 
sign in the definition of vf{t) this means that Wi is reset from +oo to — oo and the 



other variables are incremented by Js- Hence w(t) is a solution and (2.1) has time 
reversal symmetry. 

We also show that the splay state is transformed to itself by the time reversal. A 
splay state is a solution v{t) characterized by the following properties 

v,{t + T)^v,+i{t), v,{t + NT) = v,{t) j = l,...,N. (5.1) 

Note that if w{t) is the time reversal of v(t) then Wj{t) ~ —VN-j+ii—t), j ~ 1, . . . , N. 
We now make the following computation: 

Wj (t + T) = -VN-j+ii-t ~ T) 

= -VN-ji-t-T + T) (5.2) 

It follows that w{t) is also a splay state. Moreover, by choosing the phase, we can 
set wi(0) = 0, which implies that ■i;i(0) = W7v(0), or vi{0) = widN - 1)T). Therefore 
w(t) must be a phase shifted version of v(t). 
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We now use the following well known result [21] 
Theorem 1. Let 



X e 7^ 



N 



(5.3) 



be an ordinary differential equation and R a matrix. 



Suppose that (|5.3|) has a time re- 



versal symmetry defined as follows: ifx{t) is a solution of (5.3) theny(t) = —Rx{~t) 
is also a solution. Suppose also that (5.3 1 has a periodic solution Xo(^) such that 
— i?xo(— = Xo(t + T), for some T. Then all the Floquet multipliers of yL^it) are on 
the unit circle. 

It follows from Theorem [T] that the splay phase solution has all its Floquet mul- 
tipliers on the unit circle (as shown in Fig. 5.1). In particular, in Fig. 5.1 we report 



the Floquet multipliers for two different shape of PSP, but maintaining the the same 
coupling weight G, we observe that the multipliers which were inside the unit circle 
attain modulus one by passing continuously from step to (5-pulses. 




-1.0 -0.5 0.0 0.5 1.0 



Re(|x) 



Fig. 5.1. Floquet multipliers for splay state with different PSPs: namely, blue stars refer to 
step functions with J = 10, and red circles to 5-functions. The coupling weight is the same in the 
two cases, G = 180 ms. 



6. Continuous Family of Periodic Solutions. We want to show that the 
— 3 directions of neutral stability for the splay state are not only local but also 
global. We have verified this issue numerically, by perturbing randomly the splay state 
X and by following the system dynamics, with the aid of the general event-driven map 
discussed in §3A, until its convergence to some stationary state. In particular, the 
initial conditions for these simulations have been generated as follows 

x = x + aj\f , (6.1) 

where x identifies the splay state. A/" is a iV-dimensional random vector whose com- 
ponents are (5-correlated with zero average and Gaussian distributed with unitary 
standard deviation, and the noise amplitude is cr = 0.1. By following the time evolu- 
tion for a sufficiently long time span (typically, of order of 100 • N spikes), we always 
observe that these initial conditions converge to periodic orbits or to the quiescent 
state X = {—1, —1}- Tliis has been verified for system size up to A = 1, 000 and 
by considering up to 10, 000 different initial conditions for each N. 

Furthermore we observe that the final state is an orbit with periodicity x = 



if > 4 and periodicity x = 2 if TV = 4 (Fig. 6.1). For N = 3 the final state is 
always the splay state. These solutions are characterized by neurons firing periodically 
with the same period, but with time intervals among successive firing which are not 
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constant, like for the splay state. Please notice that, in the event-driven map context, 
the splay state amounts to a fixed point of the dynamics. 

All the the periodic orbits we found lie on the {N — 3)-manifold associated to 
the neutrally stable directions of the splay state in the event-driven map formulation. 



that can be obtained by eq. (4.20). We can affirm this, since on one hand we have 
verified that by perturbing the splay state along the stable directions we end up to 
the splay state itself (as shown in Fig. |6.1[ a)), while while by perturbing along the 
neutrally stable directions we always end up in one of these many periodic orbits (as 
shown in Fig. |6.1[ b)). On the other hand, by perturbing one of these orbits along the 
stable directions of the splay state the perturbed system converges to the same orbit 
(see Fig. |6.1[ c)), while by perturbing along the neutrally stable directions the system 
ends up on a different periodic orbit (see Fig. |6.1[ d)). Therefore these periodic orbits 
are also neutrally stable and share the same neutrally stable manifold of the splay 
state. 

The existence of this manifold made of a continuous family of periodic solutions 
has been previously reported for Josephson arrays [27kd3\ and Watanabe and Strogatz 
discussed the generality of this issue, reporting a "heuristic" argument to support the 
existence of this manifold for generic fully coupled oscillator networks with sinusoidal 
coupling [29]. 

As a last point we have evaluated for the splay state and several periodic orbits 
(namely Nt = 10, 000) the single neuron firing rate i^. The distribution of these rates 
is reported in Fig. |6.2[ revealing that the splay state is characterized by the minimal 
firing rate with respect to the ones found for the associated family of periodic orbits. 

7. Conclusions. In this paper we showed analytically that finite-size all-to-all 
pulse-coupled excitatory networks of excitable neurons admit marginally stable per- 
sistent splay states. We obtained analytical information about the stable firing rates 
of these sustained activities. Since the firing rate of persistent states is an electrophys- 
iologically measurable quantity in a working memory tasks, these results can provide 
insights for working memory models. We further obtained results on the splay state 
stability that can help in choosing the correct parameters required for biologically 
relevant working memory models. Our results also give an analytical basis to previ- 
ous observations in models that the stable sustained neural activity is asynchronous 

[201 [g. 

We developed event-driven map methods to analyze the network dynamics and 
found an analytical expression for the Floquet spectra associated to the splay state 
for step pulses and ^-pulses. In the case of M overlapping synaptic step pulses our 
analysis has revealed that for a correct treatment of the linear stability analysis, the 
evolution of M additional variables, corresponding to the last M firing events, should 
be taken in account. 

Our analysis, extending previous results for systems with sinusoidal coupling [29], 
revealed that the splay state is marginally stable for finite size networks with N — 2 
neutral directions, which reduce to ~ 3 in the event-driven map formulation. We 
also reported a rigorous proof for non overlapping step pulses. We further identified a 
continuous family of periodic solutions surrounding the splay state. Their peculiarity 
is that these periodic states have exactly the same neutral stability directions as the 
splay state. 

Our results leave several open questions, in particular we need to proof that 
at least one of the splay states is Lyapunov stable, when they exist. It would be 
also of interest to extend the rigorous results reported in §C to overlapping PSPs. 
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number of spikes number of spikes 

Fig. 6.1. Examples of trajectories (red lines) emerging from the perturbation of the splay state 
(black lines) or of a periodic state (blue lines). Only the voltage variable x\ is reported here as a 
function of the index labeling the sequence of successive firings. Perturbation of the splay state: (a) 
along the directions of stability, the system converges to the splay state; (b ) along the directions of 
neutral stability, the system is set in a periodic state. Perturbation of a periodic state: (c) along the 
directions of stability, the system converges to the periodic orbit; (d) along the directions of neutral 
stability, the system is set in a new periodic orbit. The system parameters are N = 5, J = 15, and 
Ts = 3 ms, and cr = 0.2. 



Furthermore, since the stable persistently active solutions of our network have a spe- 
cific spiking structure, splay or families of periodic solutions, it would be interesting 
to identify the structure of the unstable states that form the separatices between 
this sustained activity and the ground state. Finally we should explain why all the 
marginally stable states are periodic. 
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Appendix A. Explicit solution of the splay state firing rate for small 
network sizes. In this Appendix we will show how it is possible to obtain explicitly 
the firing rate v of the splay state, for iV = 2, 3, 4, and for M = (namely for J < Ji). 
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Fig. 6.2. Probability distribution of the single neuron firing rate v for the periodic solutions and 
the splay state (vertical black line). The Nt initial condition are generated by randomly perturbing 
the splay state along the directions of neutral stability, the perturbations were gaussian distributed 
with zero average. The model parameter are the same as in the previous figure, and Nt = 10,000. 



A.l. Step pulses. Eq. (3.15) can be rewritten in the following way 

(„)_ -(1-7) + (1+7K["^ 



y. 



(I + 7) -(I -7)^1"^ 



(A.l) 



where we have made use of the variable 

7 = exp(-2Ti/r) ; (A.2) 
Since in the present case < Tg < T the values of 7 are bounded between and 1. 



By employing (3.14), (A.l), and (3.18), the coefficients of the event-driven map 



(3.19) can be rewritten as 



ao = -(l-7) + (l+7)(^-l)/3im) 
fli = (l + 7) + (l-7)/3i(T,) 

a2 = (l + 7)-(l-7)(J-l)^i(r.) 

a3 = -(l-7)-(l + 7)/3im) 



(A.3) 



The firing rate can be obtained in an explicit form by inverting (A.2), namely 

1 



1 



1 



NT ATT,- f ln(7(iV, J,T,)) 



(A.4) 



Once fixed the network parameters, an admissible solutions for 7 € [0; 1] amounts 
finding a splay state solution with a frequency given by (A.4 1. 



Given an admissible 7 value, the membrane potentials corresponding to the splay 



state can be found by iterating the map (3.19) starting from the boundary condition 
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x{N) = —oo corresponding to the reset value, namely 

XN-2 [aaaz + a{) / {az{ai + a2)) 



\ 



X2 -(floaa + a|)/(a3(ai + 02)) 

xi -02/03 

Xq J \ 00 ) 



(A.5) 



We can finally determine v analytically for N — 2, 3, 4. 
• In the case of a couple of neurons, N = 2, we should impose ii = xat-i and 
thus we have ai + 02 = 0. Solving this equation for 7 we obtain: 



^ (J-2)/3i(r,)-2 
^ (J-2)/3i(T,) + 2 ^ ■''^ 

in this case we have a unique stable branch of solutions, as shown in Fig. |3.2| 
Furthermore, the minimal reachable frequency is zero and it is achieved for 
7 = 0, when J and Ts satisfy the equation (J — 2)l3i[Ts) = 2. 
For = 3 we have xi — xn-i, using the values in (A.5) we obtain: 



0903 + + a2 + 0102 = 



(A.7) 



and then, we can reorder this equation as a second order equation for 7: 



[( J- 2)/3i(T,) + 21272 - 2[( J2 - 2J + 2)/?2(T,) - 2]7+ [( J- 2)/3i(T,) - 2]^ 

this equation admits the following two solutions 
71,2 = {[(J2-2J + 2)/32(T,)-2]± 

±V[(J2 -2J + 2)Pl{T,) - 2]2 - [(J - 2)2/32(r,) - 4]2}. 

•{[(J-2)/32(T,) + 2]2}-i 



(A.8) 



(A.9) 



7i (resp. 72) is associated to the upper stable (resp. lower unstable) branch 

reported in Fig. |3.2| In this case the upper branch is bounded away from 

the zero frequency, and the minimal frequency is attained for 71 =72, when 

J and Ts satisfy (J^ — 3 J + 3)/32(Ts) = 3. The zero frequency is instead 

reachable on the lower branch for 7 = as shown in Fig. |3.2| 

If A^ = 4 then X2 = X]s[^2 and the coefficients should satisfy the following 

equation 



2aoa3 + flj 







(A.IO) 



Similarly to the case A^ = 3 we obtain a quadratic equation for the parameter 
7, namely 

[(J-2)/32(T,)+2]272_2[j2/32(j.^)]^^[(j_2)/3i(T,)-2]2=0 , (A.ll) 

also in this case we have 2 branches of solutions for the splay state frequencies 
parametrized by 71 and 72 



^ j2/32(r,) ± ^[J^Pl{Ts)f - [{J - 2YPl{Ts) 
[(J-2)/32(T,)+2p 



(A.12) 
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Also in this case the zero frequency is attainable on the unstable branch 
for 7 = and the merging of stable and unstable branch occurs at a finite 
frequency corresponding to a value of J which is solution of ( — 2 J + 
2)PKTs) = 2. 



A. 2. (5-pulses. We can rewrite the coefficients (3.24) of the map (3.19) for the 
case of (5-pulses combining (3.22), (A.l) and ( |3.18| ) as follows 



ao = -{1 --f) + {! + -/) Js 
ai = (1 + 7) 

02 = (1 + 7) - (1-7)^5 
a-i = -(1 - 7) 



(A.13) 



where 7 is given by the expression ( A.2[ ) with Ti = T. The firing rate for the splay 
state can be obtained from the following expression 

(A.14) 



^ln(7(iV,G,T)) 



Let us now discuss of the existence of the splay state for for = 2, 3, 4: 



In the case of a couple of neurons, 
obtains: 



N — 2, solving this equation for 7 one 



7 



Js 



(A.15) 



like in the step pulses case one has only one branch and the splay state exists 
for Js > 2 and the period diverges to infinite at Js — 2. 
If TV = 3 



72 



11.2 



-2±2^/jf^ 
{Js + 2)2 



(A.16) 



now two branches are present and similarly to the step pulses the upper 
branch (corresponding to 71) is stable while the other one is unstable. The 
branches exist for Js > a/S and they merge exactly for this coupling value. 

For A^ = 4 



7l,2 



Jj ± 2^ 
{Js + 2)2 



(A.17) 



also in this case the two branches are present above a certain critical coupling 
given by Js = \/2. 

Appendix B. Analytic expression for the splay state in the infinite size 

limit. In the limit of A cx) it is possible to derive an analytic expression for the 
membrane potentials associated to the splay state both for step and i5-pulses. In such 
a limit the mean input current / can be assumed to be constant, and it can be easily 



obtained from (2.4), giving / 



1. Thus we can rewrite (2.1) as follows 



dv 



222 

TT T U 



We can then integrate equation (B.l) between the reset value v 
generic time ti. 



(B.l) 



-00 and a 
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the integration gives: 



v{ti) = — 7rri/tan(- 



NT 



U) 



(B.3) 



where for the splay state v = 1/ (NT) . If we identify U with the spike time of neuron i 
in the network, we will have that the splay state solution for the membrane potential 
of neuron i is Xi = —v{ti), please notice that potential Xi are ordered from the largest 
to the smallest. Furthermore, since the spike times are equally spaced for the splay 
solution as t, 



iT, i = 0,1, ...N we can rewrite (B.3) as 



TITV 

tan(7r^; 



N- 



^ m = 



tan(7r^) 



(B.4) 



where < ^ < 1 is a continuous spatial variable. As shown in Fig. |B.H the expres- 
sion obtained in the continuous limit compare reasonably well with the numerically 
estimated finite size solutions already for TV = 16. 




Fig. B.l. Membrane potential values as a function of ^ = i/N for a splay state. The symbols 
refer to N = 16, while the solid line to the continuous limit approximation. The data have been 
obtained for J = 15, Ts = 1 ms and t = 20 ms. 



Appendix C. Marginally stable directions of the splay states. In this 
appendix, we will analyze the stability of a splay state in the case of non-overlapping 



step pulses. To perform this analysis, let us rewrite the QIF model (2.1) as follows 



_d(h 
dt 



I{t) + (/(t) - 2) cos 6ii 



1, 



(C.l) 



where we have performed the transformation of variable di = 2tan~^{vi). Therefore, 
the membrane potential is now represented by a a phase variable 9i e [— 7r;7r], the 
spike is emitted (and transmitted instantaneously to all the neurons in the network) 
whenever Oi reaches the threshold tt and then tt it is reset to — tt. The model in the 
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formulation (C.I) is termed O-newcon, we will apply the Watanabe and Strogatz j29) 



approach to this model to derive the Floquet spectrum for the splay state solution. 



In order, to stress the peculiar PSPs we are considering, we rewrite (C.ll as 
follows 

^^{jm~2)cos{9,) + jm j = l,...N (C.2) 

with 4>{9) being the characteristic function of the interval [^t:,0off]- The emission 
of a spikes occurs whenever the neuron min{(?i} — —tt this amounts to an increase 
by one in the value of the function (j){9). Furthermore, when the pulse expires after a 
time Tg the value of (f){9) will be decreased by one. By assuming that no neuron will 
fire while the synapse is on (no overlapping PSPs), the PT will occur for a specific 
value of the phase variable namely mm{9i} — 9off, for a value Oqff, which can be 
determined as outlined in §3A. 

Let us now recall the approach devised by Watanabe and Strogatz [29] to show 
that each trajectory representing the dynamics of a system of N identical phase 
oscillators, whose evolution is ruled by ODEs of the form 

'^ = f{e)+g{e)cos{9,) + hi9)sm{0,) j = l,...iV , (C.3) 

is confined to a three dimensional subspace. The only requirement is that the functions 
/, g and h do not depend on the index j of the considered oscillator. In other words /, 
g and h are collective variables determined by the network state. Clearly our equation 
(|C.2|) satisfies this condition. 



Watanabe and Strogatz introduce a transformation : — > R^^^ from vari- 
ables {9j} to variables X = (F, 0, 5', {^pj}) defined implicitly by the set of equations 

i^(0„F,e,vI/,V'j) = O, j = l,...N ; (C.4) 

where 



F = tan ( ^-{9, - 6)] - J^^t^n - ) . (C.5) 



Furthermore, they prove that an arbitrary solution of (C.3 1 can be generated by the 
transformation Qx from a set of parameters {ipj}i which remain constant in time, 
whenever the three collective variables F, 9, ^' satisfy the following equations 

f = -(1 - F2)(gsine - /icosG) 

Fe = F/-gcose-/isine (C.6) 



F* = -r2(5COse + /isine) 
and obviously the other variables satisfy 

^,=0 Vj = l,...,iV . (C.7) 
We prove the following proposition: 

Proposition C.l. Let us assume that (C.2) admits a splay state solution and 
that this solution is Lyapunov stable. Then at least N — 2 Floquet multipliers will lie 
on the unit circle. 
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Let us now recall the definition of the Floquet multipliers j7j for a generic ODE 
of the form 



e R 



N 



(C.8) 



admitting a periodic solution 0''{t) with period Tp. 

The associated variational linear equation in the tangent space has the form: 



60 = DF,{9'it))Se, 



e n 



N 



(C.9) 



Equation (C.9) has (possibly complex) eigensolutions <I>(t) = e^-^~^'^'^'>*r]{t), with ri{t) 
periodic of period Tp, termed Floquet vectors. The complex numbers fJ.{Tp) = 
^{x+iuj)Tp g^j.g i-j^g Floquet multipliers. They determine the stability of the periodic 
solution. 



Proof of C.l We will prove that there exists an — 2 dimensional subspace 
of solutions of the variational equation associated to (C.2| consisting of solutions 
that do not converge to the vector &s t ^ co (except for the solution itself). 
This, combined with Lyapunov stability, implies that there must be iV — 2 Floquet 
multipliers on the unit circle. 

Let 6'g = {6*0,1; • • • 7 ^0 Af} be a choice of initial conditions corresponding to a splay 
state of period Tp. For simplicity and without any loss of generality, we can assume 
that the phases are ordered, i.e. i > ^0,2 > • ■ • , %.n^ ^'^^ ^^^^ ^0.1 close to tt, i.e. 
the first neuron is just about to fire. 

ition X''(t 

n = {0, 



Let us consider a solution X''(t) of (C.6) and (C.7) with initial condition 

TT TT 

2' 2' 



^0,1' ■ 



} 



(C.IO) 



where it is evident that 6'"(t) = T:,{X.''{t)) since Tr(X"(0)) = 6''*(0). Furthermore, we 
perturb the initial condition with a perturbation of the form 



A^ = (Ai^i,...,AV'jv-2,0,0) 



(C.ll) 



in the following way 
X1^(0)=X§ + AV 



(C.12) 



and we obtain the perturbed solutions X^^(i) at time t by integrating (C.6) and 
(C.7 1, while the the corresponding solution of ( C.3[ ) is given by 0^^^{t) — Tx{X.''^^{t)). 
Let us denote the value of the perturbed orbit at integer multiples k of the period 



Tp as follows 



= inkTp),e{kTp),^j;{kTp),9l + A^, 



and 9 



Atl:,k 



(C.13) 



We will show that there exists a real positive constant L such that, for every value 



of fc, 



>L\\Ai^\\. 



(C.14) 
-Xgll « l,we 

can approximate the evolution of the perturbed orbit in proximity of the unperturbed 
one, with the corresponding linearized dynamics, namely 



By assuming that the perturbation is sufficiently small, i.e. | |X 



A4:,k' 



T,(Xl^,)-T,(Xg)«OT.(X^)(X 



AipM 



xg) 



(C.15) 
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In order to write the Jacobian DTx(X.q), we need to estimate the foUowing derivatives, 
which can be obtained by imphcit differentiation of (C.4| 



d9^ 

de 



= 1, 



9* 



= -1, 



89^ 



^-6 



k- 



(C.16) 



where 5jk is the Kronecker delta. 

Let V° S TZ^ be a vector with an unitary norm spanning a — 2 dimensional 
subspace and let assume that Ai/i — ctV" with < cr < 1. We will prove that 



1 



i?r,(x^)(x 



Xg)|| >L>0 



(C.17) 



for some real constant L independent of a and k. 



By employing (C.16) the following expression can be derived 



/ Zi(fc)/a + yO 
Z2{k)/a 



-OT,(xg)(X: 



Aip.k 



ZN-2{k)/(J + Vji_^ 

ZN-i{k)/(j 
ZNik)/a 



(C.18) 



where for brevity and clarity we set Zj{k) = — z^. — cosO^v^. once redefined = 
T{kTp), Wk = Q{kTp) — |, and Zk — ^'(feTp) — ^. It is clear, due to their definition, 
that the components of the vector Z(fc) = {Zj[k)} are not linearly independent and 
in particular that they span a 2-dimensional subspace. 

As a first step, the validity of the following inequality, V/c and for any sufficiently 
small a, is discussed 



|Z,(fc)|/a ^\{wk-Zk- cose'vk)\/a > L for j ^ N or j ^ N - 1 



(C.19) 



We consider two possible cases. In the first case, the inequality (C.19 1 holds 



therefore (C.17 1 is satisfied since the length of any vector is bigger than the absolute 



value of one of its components, thus implying that the modulus of the l.h.s. of (C.18 1 
would be greater than L for any k value. 



In the second case, we assume that (C.19) does not hold uniformly in k for j — N, 
and J = iV — 1, in other words the components |Zjv_i(fc)|/(T and \ZN{k)\/cr should con- 
verge to for fc — >• cx) and a oo. Furthermore, since for > 3 cos0^ ^ cos6ff_j^, 
each component Zj with j — 1, . . . , — 2 can be written as a linear combination 
of Zjv_i and Z^. This implies that each element |Zjv-i(fc)|/o' remains arbitrarily 
small Vj even for arbitrarily large (resp. small) k (resp. cr). Now each component 
in the r.h.s. of ( |C.18[ ) will have the form Zj/a + V° for j ^ 1, . . . , N - 2, where the 
first quantity is arbitrarily small, but by construction the vector V*^ has an unitary 



modulus, thus also in this second case (C.17) is satisfied for any k. 



From the previous results it follows that the vector function 



V(<) = — 0A*(i)U=o 

da 



(C.20) 



is a solution of the variational equation (C.9 1 which does not converge to as t — )■ oo 



Since (C.9 1 is a system of linear equations, a vector space of initial conditions gives 
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rise to a vector space of solutions. Since V'^ spans a — 2 dimensional vector space, 
which we denote by LV, our construction give a iV — 2 dimensional vector space of 
solutions of (C.9|, which we denote by £V. 



As mentioned above the Floquet vectors are solutions of ( C.9 ) of the form ^{t)rj{t), 
with 'q{t) periodic of period Tp and ii{Tp) the corresponding Floquet multipliers. Since 
we assumed that the examined periodic orbit (i.e. the splay state) is Lyapunov stable, 
the multipliers n{Tp) must be either on the unit circle or inside the unit circle. Without 
loss of generality, let us assume that at least two multipliers are inside the unit circle, 
otherwise the theorem would be automatically true. 

Let denote by LW the vector space spanned by the initial conditions of the two 
Floquet eigenvectors associated to the two multipliers which lie inside the unit circle 
and let CW be the corresponding vector subspace of solutions (spanned by the two 
Floquet eigenvectors). Since all non-zero solutions in £W converge to as t — >■ cxi it 
follows that the intersection of LW and LV consists of the zero vector. Therefore, 
we can formally decompose any of the remaining N — 2 Floquet vectors at initial 
time i = in two vectors, namely 7^(0) = Wi(0) + Vi(0) where Wi(0) € LW and 
Vi(0) e LV. By linearity, if Wi(t) and Vi(t) are the solutions of (C.9 1 with initial 
conditions Wi(0) and Vi(0), it follows that T]{t) = Wi(t) + Vi(i). If Vi(0) ^ 
then ri{t) ^ CW, moreover ri{t) does not converge to as t ^ oo since Wi(t) does, 
while Vi{t) does not. Therefore the corresponding Floquet multiplier can be only on 
the unit circle, due to our previous assumptions. Finally we have demonstrated that 
N — 2 Floquet multipliers are on the unit circle and 2 are inside the unit circle. 
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